package com.googlecode.curvedmetromaps.core.geometry;

public class BSplineFitBasis {
	// The number of control points and degree for the curve.
	private int			mQuantity, mDegree;

	// The storage for knots and basis evaluation.
	private double[]	mValue;			// mValue[0..degree]
	private double[]	mKnot;				// mKnot[2*degree]

	// Construction and destruction. This class is only for open uniform
	// B-spline basis functions. The input is the number of control points
	// for a B-spline curve using this basis and the degree of that curve.
	public BSplineFitBasis(int quantity, int degree) {
		assert 1 <= degree && degree < quantity : "Invalid inputs.";
		mQuantity = quantity;
		mDegree = degree;

		mValue = new double[degree + 1];
		mKnot = new double[2 * degree];
	}

	// / Getter for mQuantity
	public int GetQuantity() {
		return mQuantity;
	}

	// / Getter for mDegree
	public int GetDegree() {
		return mDegree;
	}

	double GetValue(int i) {
		assert 0 <= i && i <= mDegree : "Invalid index";
		return mValue[i];
	}

	void Compute(double t, int[] minMax) {
		int imin = minMax[0], imax = minMax[1];
		assert 0.0 <= t && t <= 1.0 : "Invalid input.";

		// Use scaled time and scaled knots so that 1/(Q-D) does not need to
		// be explicitly stored by the class object. Determine the extreme
		// indices affected by local control.
		double QmD = (double) (mQuantity - mDegree);
		double tValue;
		if (t <= 0.0) {
			tValue = 0.0;
			imin = 0;
			imax = mDegree;
		} else if (t >= 1.0) {
			tValue = QmD;
			imax = mQuantity - 1;
			imin = imax - mDegree;
		} else {
			tValue = QmD * t;
			imin = (int) tValue;
			imax = imin + mDegree;
		}

		// Precompute the knots.
		for (int i0 = 0, i1 = imax + 1 - mDegree; i0 < 2 * mDegree; ++i0, ++i1) {
			if (i1 <= mDegree) {
				mKnot[i0] = 0.0;
			} else if (i1 >= mQuantity) {
				mKnot[i0] = QmD;
			} else {
				mKnot[i0] = (double) (i1 - mDegree);
			}
		}

		// Initialize the basis function evaluation table. The first degree-1
		// entries are zero, but they do not have to be set explicitly.
		mValue[mDegree] = 1.0;

		// Update the basis function evaluation table, each iteration
		// overwriting
		// the results from the previous iteration.
		for (int row = mDegree - 1; row >= 0; --row) {
			int k0 = mDegree, k1 = row;
			double knot0 = mKnot[k0], knot1 = mKnot[k1];
			double invDenom = (1.0) / (knot0 - knot1);
			double c1 = (knot0 - tValue) * invDenom, c0;
			mValue[row] = c1 * mValue[row + 1];

			for (int col = row + 1; col < mDegree; ++col) {
				c0 = (tValue - knot1) * invDenom;
				mValue[col] *= c0;

				knot0 = mKnot[++k0];
				knot1 = mKnot[++k1];
				invDenom = 1.0 / (knot0 - knot1);
				c1 = (knot0 - tValue) * invDenom;
				mValue[col] += c1 * mValue[col + 1];
			}

			c0 = (tValue - knot1) * invDenom;
			mValue[mDegree] *= c0;
		}

		minMax[0] = imin;
		minMax[1] = imax;
	}
}
